Isotope effects on the lattice parameter of cubic SiC 
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Path-integral molecular dynamics simulations in the isothermal-isobaric (NPT) ensemble have 
been carried out to study the dependence of the lattice parameter of 3C-SiC upon isotope mass. This 
computational method allows a quantitative and nonperturbative study of such anharmonic effect. 
Atomic nuclei were treated as quantum particles interacting via a tight-binding-type potential. At 
300 K, the difference Aa between lattice parameters of 3C-SiC crystals with 12 C and 13 C amounts 
to 2.1 x 10 -4 A. The effect due to Si isotopes is smaller, and amounts to 3.5 x 10~ 5 A when 
replacing 28 Si by 29 Si. Results of the PIMD simulations are interpreted in terms of a quasiharmonic 
approximation for the lattice vibrations. 
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It is well known that the lattice parameters of two 
chemically identical crystals with different isotopic com- 
position are not equal, lighter isotopes giving rise to 
larger lattice parameters. This is caused by a combina- 
tion of two factors: the dependence of atomic vibrational 
amplitudes upon atomic mass, and the anharmonicity of 
the vibrations. The isotope effect on the lattice param- 
eter is most important at low temperatures, due to the 
change of zero-point vibrational amplitude with atomic 
mass, and disappears in the high-temperature (classical) 
limit at T > 0£> (&d, Debye temperature), where vibra- 
tional amplitudes are independent of the mass. In recent 
years, it has become feasible to measure the isotopic ef- 
fect in lattice parameters of crystals with high precision. 1 
Most of the work has been performed on elemental crys- 
tals, although binary and multinary materials offer the 
attractive possibility of isotopic substitution on different 
atoms' Of these materials, SiC, with 3 stable isotopes 
of Si and 2 of C, plus about 70 different polytypes, is 
a paragon. Here, however, we confine ourselves to the 
simplest polytype: zincblende-like 3C-SiC. 

SiC has been suggested for a number of applications ex- 
ploiting many of its superlative properties, close to those 
of diamond. 3 ' 4 Some of these applications take advantage 
of its hardness, large thermal conductivity, and low ther- 
mal expansion. Isotopically modified SiC may find appli- 
cations exploiting the higher thermal conductivity 5 - and 
the dependence of its hardness on isotopic composition. 6 

Anharmonic effects in the vibrational properties of 3C- 
SiC have been studied earlier in detail, e.g., the pressure 
and temperature dependence of phonon frequencies^ 
and lifetimes^ii^ The thermal expansion of cubic SiC 
(another anharmonic effect) has been studied in detail 
both experimentally-ii and theoreticall y 12 ! 13 , the latter 
using a quasiharmonic approximation (QHA) for the lat- 
tice vibrations. 

Isotopic effects on the lattice parameters of crys- 
tals have been usually calculated by employing the 



QHA and perturbative methods based on ab-initio 
techniques & 14 ' 1 5 An alternative to perturbational ap- 
proaches in solids is the combination of the path inte- 
gral formulation (to deal with the quantum nature of 
the nuclei) with electronic structure methods. The path 
integral approach to statistical mechanics allows one to 
study finite-temperature properties of quantum many- 
body problems in a nonperturbative scheme, even in the 
presence of large anharmonicitieSfi& An advantage of its 
combination with electronic structure methods is that 
both electrons and atomic nuclei are treated quantum 
mechanically in the framework of the Born-Oppenheimer 
approximation ! 17 ' 18 ' 19 The path-integral molecular dy- 
namics (PIMD) method is based on an isomorphism be- 
tween the quantum system under consideration and a 
classical one, obtained by replacing each quantum parti- 
cle by a cyclic chain of classical particles, connected by 
harmonic springS j 16 ' 20 ' 21 

When calculating properties of crystals with isotopi- 
cally mixed composition, it is usually assumed that each 
atomic nuclei in the solid has a mass equal to the aver- 
age mass. This kind of virtual- crystal approximation has 
been used in density-functional calculations, as well as in 
atomistic simulations based on path integrals j 2 ' 21 ' 22 ' 23 ' 24 
In fact, in earlier path integral simulations of diamond 
it was found that the results obtained by using this 
approximation are indistinguishable from those derived 
from simulations in which actual isotopic mixtures were 
considered^ 

Here we extend earlier path integral calculations of 
the lattice parameter of group-IV solid o 21 ' 23 ' 24 (diamond, 
Si, Ge) to a IV-IV compound such as 3C-SiC. The 
electronic structure has been treated with an efficient 
tight-binding Hamiltonian, based on density-functional 
calculations. 25 ' 26 Simulations were performed on a 2 x 
2x2 supercell of the 3C-SiC face-centered cubic cell with 
periodic boundary conditions, including 64 atoms. For 
a given temperature and isotopic composition, a typical 
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FIG. 1: Temperature dependence of the lattice parameter 
of 3C-SiC, as obtained from (a) PIMD simulations (squares), 
(b) simulations performed by setting the Si nuclei as classical 
particles (circles), (c) simulations carried out by considering 
classical C nuclei (triangles), and (d) all nuclei are assumed to 
be classical (diamonds). Lines represent empirical fits. The 
statistical error of the simulation results is less than the sym- 
bol size. 
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FIG. 2: Temperature dependence of the isotope effect on 
the lattice parameter of 3C-SiC. Squares represent the dif- 
ference between lattice parameter of nat Si 12 C and nat Si 13 C, 
whereas circles indicate the difference Aa between 28 Si nat C 
and 29 Si nat C. Filled symbols indicate the low-temperature 
values derived from Eq. {5). Lines represent fits of the simu- 
lation results to Eq. JT}. 



run consisted of 2 x 10 4 PIMD steps for system equili- 
bration, followed by 4 x 10 6 steps for the calculation of 
ensemble average properties. Details on the actual im- 
plementation of the PIMD method to study structural 
and electronic properties of SiC were given elsewhere^ 

First of all, we quantify the influence of quantum ef- 
fects on the lattice parameter of cubic SiC. With this pur- 
pose we have performed simulations in which the atomic 
nuclei of either C or Si (or both) were considered as clas- 
sical particles. We note that in the formalism used here, 
the classical limit for a given atomic nucleus is obtained 
by taking its mass as M — > oo (in fact, we took masses 
in the order of 10 5 amu). 

In Fig.l we present the results of our molecular dynam- 
ics simulations for the various cases considered. First, we 
compare the results of the full PIMD simulations (quan- 
tum atomic nuclei, open squares) with those of classical 
simulations (diamonds). The zero-point motion induces 
an increase in the lattice parameter of 0.011 A, which 
means a relative change of 2.5 x 10 -3 . Second, we are 
interested in the lattice parameter obtained when nu- 
clei of one of the elements (Si or C) are considered as 
classical and the other as quantum particles. Thus, tri- 
angles represent results for quantum Si and classical C, 
whereas circles indicate data for quantum C and classical 
Si. By comparing these results with those found in the 
full PIMD simulations, we find for T — ► that the lattice 
parameter a increases by 3.9 x 10~ 3 A and 6.8 x 10 -3 A, 
due to having considered Si or C as quantum particles, 
respectively (assuming in each case that atomic nuclei 



of the other element are quantum particles). We note 
that these results are similar to those presented earlier 
in Ref. l26l . but here the simulations extended over runs 
10 times longer in order to improve the accuracy in the 
lattice parameter, and the present error bars are smaller 
by a factor of ~ 3. 

Results presented until now were obtained for SiC crys- 
tals with natural isotopic composition. We now turn 
to simulations of SiC crystals, in which one of the ele- 
ments is taken to be isotopically pure. Thus, we carried 
out PIMD simulations of nat Si 12 C and nat Si 13 C (where 
the superscript "nat" refers to the mass of the natural 
isotopic composition) at several temperatures. Results 
for the difference Aa between lattice parameters of both 
kinds of crystals are shown in Fig. 2 (squares). At room 
temperature, we find Aa = 2.1 x 10~ 4 A, which trans- 
lates into a relative change Aa/a = 4.8 x 10~ 5 . Some- 
thing similar has been done for 28 Si nat C and 29 Si nat C 
crystals, and the difference Aa yielded by our simula- 
tions is much lower: at 300 K we found Aa = 3.5 x 10~ 5 
A, or Aa/a = 8 x 10~ 6 . The temperature dependence 
of Aa can be fitted to a function with the shape of the 
mass derivative of a Bose-Einstein function, namely: 
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where x = T a /T, and C and T a are fit parameters. 
Dashed lines in Fig. 2 represent this function with T a = 
1200 K (upper curve) and 600 K (lower curve). Eq. (JTJ 
is basically equivalent to Eq. (3) of Ref. [26| except that 
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the adjustable parameter b has been replaced by well- 
defined physical variables more convenient for the sub- 
sequent treatment. The minor deviations of the PIMD 
points from the fitted curves are not surprising since we 
have used a single Einstein oscillator fit. They could be 
decreased by using two oscillators.— 

The results of our PIMD simulations can be further 
analyzed in terms of a QHA for the lattice vibrations. 
In such an approximation, the lattice parameter a(T) for 
a given isotopic composition at temperature T can be 
derived by minimizing the Helmholtz free energy with 
respect to the crystal volume.—^ One finds 
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where 
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Here, uj n (q) are the frequencies of the nth mode 
in the crystal, B is the bulk modulus, Goo is the 
zero-temperature lattice parameter in the limit of in- 
finite atomic mass (classical limit), and 7n(q) = 
— <91nti>„(q)/<91n V is the Griineisen parameter of mode 
ri, q. Then, at T = the difference a(0) — aoo is given by 
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Let us consider now for simplicity two isotopically-pure 
monatomic crystals with a mass difference AM. The 
difference between the corresponding lattice parameters, 
Aa, can be related to the zero-point renormalization, 
anat(O)— Ooo, for the natural crystal. This can be achieved 
through a first-order expansion for the lattice parameter 
as a function of the mass M, and taking into account 
that the frequencies w„(q) scale as 1/vM. One finds for 
the change in lattice parameter at T = 0: 



I . . . . AM 

Aa = -- (Onat(0) - floo) Jl— 



(5) 



This means that the low-temperature changes in a due to 
isotopic mass can be obtained from the zero-point renor- 
malization of the lattice parameter in the natural crystal. 

For binary compounds such as SiC, one can use a for- 
mula similar to Eq. (O to obtain the separate contribu- 
tions of each kind of atoms (i.e., Si or C). To first order, 
the contributions of both types of atoms will be additive. 
Then, from the difference a na t(0) — a^o discussed above, 
and presented in Fig. 1, we obtain in the low-temperature 
limit, using Eq. ([5]), Aa = 2.8 x 10~~ 4 A when replacing 
12 C by 13 C, and Aa = 7.0 x 1(T 5 A for substitution of 
28 Si by 29 Si. These are the values shown in Fig. 2 as 
filled squares at T = 0. 

We note that the virtual-crystal approximation has 
been employed in our simulations, i.e., for crystals with 
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FIG. 3: Lattice parameter of 3C-SiC as a function of M e g 
in a virtual-crystal approximation, in which all atoms in the 
simulation cell (Si and C) are assumed to have the same effec- 
tive mass. Open squares indicate results of PIMD simulations 
in this approximation. The dashed line is a fit to Eq. Q. 
An arrow indicates the lattice parameter derived from PIMD 
simulations with C and Si atoms having their actual mean 
isotopic masses. 



the natural isotopic composition of Si or C we have as- 
sumed that the atoms have the average mass. In this re- 
spect it is worthwhile considering the validity of assuming 
an effective mass M e ff for all the atoms in a given crystal, 
to describe changes in the lattice parameter. In view of 
Eq. ||3J), one has^ 2 . 
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where V c is the volume of the primitive cell, and (...) 
indicates an average over all branches of the Brillouin 
zone. Now we may assume a dependence of the average 
frequency on effective mass as 



(w„(q)) ~ M t 
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which, in fact, is expected when one considers M e ff as 



1 



Mr 



1 

Mc 



1 



(8) 



To check this point we have carried out PIMD simula- 
tions for SiC crystals with various effective masses for 
both Si and C. The results are displayed in Fig. 3 as 
open symbols. For comparison with the result derived 
for M e ff given by Eq. jS]) (~ 16.8 amu), we also show 
those obtained by assuming an effective mass given ei- 
ther by the average (Mc + Mgj)/2, or by the mass of Si 
or C. All these results can be fitted well to the expression 

(9) 
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where b and c are fit parameters. Such an expression 
can be expected from Eq. Q when one considers an ef- 
fective mass for both types of atoms, at temperatures 
T <C &d (for SiC, Q D ~ 1100 K). The actual lattice 
parameter of SiC yielded by the simulations above is in- 
dicated in Fig. 3 by an arrow. It coincides within error 
bars with that derived assuming M e g = 16.8 amu, as 
given in Eq. {SJ. Note that taking the effective mass as 
the average of the masses (at about 20 amu) yields a lat- 
tice parameter clearly lower than the actual one obtained 
using the separate masses of C and Si. 

The agreement between the result for M e s = 16.8 amu 
and the real crystal can be interpreted in terms of per- 
turbation theory as follows. Looking at Eq. (j4|), changes 
in the lattice parameter are mainly due to TO phonons, 
as can be derived from Fig. 9 in Ref. .13 for the appro- 
priate values of 7„ and Fig. 2 for the density of phonon 
states. The TO band in SiC is rather symmetric^ and 
according to second-order perturbation theor y 22 ' 28 mass 
fluctuations cause an increase in the high frequencies and 
a reduction of the low ones by a similar amount. All to- 
gether the effect in the lattice parameter is expected to 
be negligible, as observed in the results of the simula- 
tions. If some effect appears due to mass fluctuations, 
it has to be of third or higher order, and should be less 



than the statistical uncertainty of our results. This gives 
further support to the virtual-crystal approximation for 
calculating lattice parameters of this kind of semiconduc- 
tors. 

In summary, we have calculated the isotope effect on 
the lattice parameter of cubic SiC by PIMD simulations. 
This procedure gives a quantitative estimation of such 
effect, which amounts to Aa/a — 4.8 x 10~ 5 for replace- 
ment of 12 C by 13 C, and to Aa/a = 8 x 10~ 6 for substitu- 
tion of 29 Si for 28 Si. These results have been interpreted 
in terms of a quasi-harmonic approach for the lattice vi- 
brations. The virtual-crystal approximation is found to 
be valid in all cases considered here. Second order per- 
turbation theory, as implied by Eq. (|5|), has also been 
shown to be valid. 
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